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Abstract 

The resources required to solve the general interacting quantum A^-body problem scale exponen- 
tially with A^, making the solution of this problem very difficult when N is large. In a previous series 
of papers we develop an approach for a fully-interacting wave function with a general two-body 
interaction which tames the A^-scaling by developing a perturbation series that is order-by-order 
invariant under a point group isomorphic with Sm . Group theory and graphical techniques are 
then used to solve for the wave function exactly and analytically at each order. Recently this 
formalism has been used to obtain the first-order, fully-interacting wave function for a system of 
harmonically-confined bosons interacting harmonically. In this paper, we report the first applica- 
tion of this A^-body wave function to a system of N fully-interacting bosons in three dimensions. We 
determine the density profile for a confined system of harmonically-interacting bosons. Choosing 
this simple interaction is not necessary or even advantageous for our method, however this choice 
allows a direct comparison of our exact results through first order with exact results obtained in 
an independent solution. Our density profile through first-order in three dimensions is indistin- 
guishable from the first-order exact result obtained independently and shows strong convergence 
to the exact result to all orders. 
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I. INTRODUCTION 



The interacting quantum A^-body problem becomes particularly challenging when N is 
large. Unlike the corresponding classical problem where the resources required to solve the 
problem scale as a polynomial in , the resources needed to solve the quantum A^-body 
problem scale exponentially with A^, frequently doubling for every particle added. Q] 
When interparticle interactions are weak, the mean-field approximation may be used to 
avoid this exponential scaling with A^ . When interactions are larger, phenomenological 
models are often used. Typically the phenomenological models are only valid for a range 
of interaction strengths, masses, . . .etc., so more general first-principles approaches which 
tackle the quantum mechanical exponential A^-scaling issue head on are needed. Such 
first-principles methods for confined, A^-particle quantum systems include coupled cluster 



methods (CCM) 



^, the method of correlated basis functions (CBF)0, H], and density 



functiona 
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theory 

n 
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i4M,m- 

We take a different approach and develop a non-numerical method which uses group 
theory and graphical techniques to tackle the A^-scaling problem (see Ref. [iT] and the refer- 



ences therein). This method circumvents the severity of the A^-scaling problem by using a 
perturbation expansion about a maximally symmetric structure in large dimensions which 
has a point group isomorphic to Sjy . Group theory is then used to separate the A^ scaling 
problem away from the interaction dynamics, allowing the A^ scaling to be treated as a 
straight mathematical issue. The perturbation expansion of the Hamiltonian is order-by- 
order invariant under the point group, yielding a problem at each perturbation order 
that can be solved, essentially exactly and analytically, using group theory and graphical 
techniques. As part of this solution, small finite, S'Ar-invariant basis sets are used that are 
complete at each order, and as A^ increases, the group theory and graphical techniques "hold 
their own" with the result that the number of elements of this basis do not grow with A^ . 
(The basis elements must remain invariant under the A^! operations of the Sn group which 
puts increasing restrictions on the growth of the set as A^ gets larger.) The completeness of 



17. Since the elements of this basis, called 



this basis at each order was established in Ref. 
binary invariants, are themselves invariant under the maximal point group symmetry, the 
invariance of the Hamiltonian at each order is naturally ensured by expressing it in terms of 
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this relatively small basis (seven elements at lowest order in the wave function, twenty-three 
elements at next order). 

Applying this approach at lowest order, we have previously derived beyond-mean-field 



energies 



18 



sity profiles 



19], frequencies [18^, normal mode coordinates 201]. wave functions 201] and den- 
211 ] for general isotropic, interacting confined quantum systems. More recently. 



17 and 



22, we have extended this anal- 



in a major development of the method presented Refs. 
ysis to first order in the wave function, which required the development of new techniques to 
handle, exactly and analytically, the complexity of coupling the many normal coordinates of 
the lowest-order solution. This basic approach, developed at first order, is general enough to 
be extended to higher orders in the perturbation expansion. In Ref. ]l7), the general theory 



of Ref. 



22 



3, which derives the wave function through first order for an arbitrary isotropic 
system, was tested on the exactly soluble system of harmonically-interacting particles 
under harmonic confinement. When this wave function was compared to the exact analytic 
wave function obtained in an independent solution and then expanded analytically through 
first order, exact agreement was found, confirming this general theory for a fully interacting 
iV-body system in three dimensions. [1 71] 



22 



exact 



In the present paper, we test this general, fully interacting wave function of Ref. 
through first order, by deriving a property - namely the density profile which is directly 
observable in the laboratory for a Bose-Einsten condensate. We test it on the same exactly 
soluble system of harmonically- interacting particles under harmonic confinement. Our 
density profile through first-order evaluated at D = 3 is indistinguishable from the D = 3 
first-order result from the independent solution and shows strong convergence to the exact 
D = 3 result to all orders. 

This derivation of the density profile through first order builds upon previous isotropic 
work. The isotropic, lowest-order ground state wave function was derived in Ref. l20| and 



the corresponding lowest-order ground state density profile was derived from this in Ref, 



The isotropic, first- order ground state wave function was derived in Ref. 



in Ref. 



22 



and checked 



171 . This work is reviewed in Section [TTl while Appendix [A] contains the briefest 



review of binary invariants. The derivation of the density profile through first order from 
the general wave function through first order is found in Section UTTl Section HVl presents the 
application and results. Appendix [B] presents the exact independent solution for the density 
profile through first order from the full D-dimensional density profile for this system, which 
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is also derived independently in Appendix [B] from the ful 
the harmonically interacting system discussed in Ref. Il7 . 
conclusions section. 



D-dimensional wave function for 
Section |V] is the summary and 



II. THE GROUND-STATE WAVE FUNCTION 



A. Lowest-Order Wave Function 



The zeroth-order Hamiltonian is that of a mult i- dimensional harmonic oscillator (see 



Eq. (20) of Ref. |2J, or Eq. (4) of Ref. llTl ). Thus upon transformation to the normal modes 



of the system, the wave function is a product of P = A^(A^ + l)/2 one-dimensional harmonic 



oscillator wave functions 



20 



221. 



(1) 



u=l 



where Ui, is the frequency of normal mode ql , and rii, is the oscillator quantum number, 
< n,y < (yo, which counts the number of quanta in normal mode u. 



B. Normal Modes and Symmetry Coordinates 

The transformation to normal modes would appear to be a formidable proposition since if 

is in the millions, the number of normal modes P is of the order 10^^ or larger. However, 
the D ^ oo structure is maximally symmetric (see Ref. Il9l) . and it is this maximal point 
group symmetry which allows the normal modes to be derived (see Ref. I20i ). 

Under this point group symmetry, the normal modes transform under the [A^] , [A^ — 
1, 1] , and [A^ — 2, 2] irreducible representations. Two normal modes transform under two 
one-dimensional [A^] irreducible representations (irreps.) of Sn , 2(A^ — 1) normal modes 
transform under two (A^ — l)-dimensional [A^ — 1, 1] irreps., and A^(A^ — 3)/2 normal modes 
transform under an A^(A^ — 3)/2 dimensional [A^ — 2, 2] irrep. . 

The normal modes may be written as 

q'l = cl (cos ei [5?,]^ + sin 6% \S^,]^) (2) 

for the a = [N] and [A^ — 1, 1] sectors, where the ± denote the two normal mode vectors 
for each a , and 
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/[7V-2, 2] ^ [N-2, 2] 2] 



(3) 



where the symmetry coordinates [S'^/]^ (defined in Ref. l20l ) may be written 

N 



[iV-l,l]n 
r' \ 



fc=l 



N{N -I) 



N l-l 
• =2 A:=l 



ir' 



^,(, + l)(iV-2) 
where 1 < i < — 1 , and 



i+l J ' 

i l-l 



N 



1=2 k=l k=l l=k+l . 



— t 



N 



k=l 



1=1+2 



,^IN-2,2U 



1 



v/^(« + l)(j-3)(j-2) 



Y Y%' + Y Y -^)Y%- 

J'=2 k=l k=l j'=k+l k=l 

Y^'ui+i) + Y ) + - 3)7(i+i)j 

vfc=l j'=i+2 / 



where 1 < i < j — 3 and i + 3 < j < . 



(4) 



All of the normal modes belonging to a normal mode vector transforming irreducibly 
under Sn have the same collective motion frequency, so that instead of A^(A^ + l)/2 possible 
distinct frequencies we only have five distinct frequencies. Thus writing the wave function 
in terms of the abbreviated irrep. labels = [A^] , 1 = [A^ — 1, 1] , and 2 = [A^ — 2, 2] , 
we have 

$o(qO= n n '^"(M>0(v^[ik) , (5) 

/.={0±, 1±, 2} ?=1 

where 0n(^,5) {y/^ Wfils) ^ one-dimensional harmonic-oscillator wave function of frequency 
ujf^ and n(/i, ^) is the oscillator quantum number, < n(/i, ^) < oo , which counts the number 
of quanta in each normal mode. The index /i labels the manifold of normal modes with the 
same frequency ci)^ while degeneracy of the fith normal mode is denoted rf^ = 1 , A^ — 1 or 
N{N - 3)/2 for = 0± , 1± or 2, respectively. 
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C. Lowest-Order Ground-State Wave Function 



The lowest-order wave function po{c^) for the ground state is given by Eq. ([5]) with all 
of the rii, set equal to zero, i.e. 

p 

^o(q') = n '^o (v^ ^0 ' (6) 

i/=i 

where (pQ is the wave function of a single harmonic-oscillator, 

(po{V^ql) = (^^)^ exp (-^u^q'Jj . (7) 



D. First-Order Wave Function 

Previous applications of dimensional perturbation theory went to very high order in 
the asymptotic 1/D expansion. For systems with a large number of degrees of freedom, 
the derivation of these high order terms can be computationally prohibitive and subject 



to numerical difficulties. In Ref. 



23l Dunn et. al. present an algorithm by which these 



corrections may be derived exactly using tensor algebra. Using this formalism, the first- 

22l . Writing the wave 



order correction to the lowest-order wave function is derived in Ref. 
function through first-order as 

$i(q') = (l + 5U)$o(q'), (8) 
then A satisfies the commutator equation 

[A,i?o]$o = i?i$o. (9) 

To solve this equation, it is helpful to note that since $o(ct') is a Gaussian function, the 
derivatives in Hi and Hq written in normal coordinates "bring down" normal coordinates 
from the exponent so that Hi effectively becomes a 3rd-order polynomial of only odd powers 
in q'. 

1. Evaluation of Derivatives in Hi $0 

We evaluate the derivatives implicit in Eq. ([9]), noting the Gaussian form of $o(q')5 
reduce the operator equation to a polynomial equation in , 

d,'M^) = -^uQ'.M^) (10) 



%$o(q') 
Therefore with the substitutions 



(11) 



(12) 
(13) 



the action of Hi on $o(q') becomes equivalent to the action of a 3rd-order polynomial (^Hij^^ 
on $o(q'): 

i?i$o(q') = (i?i),g$o(q') (14) 

where 



1 



(1) 



3 



1 r 



1^1,1^2,1^3 



t^^2 + - 



'^l,l^2,'^2 



(1) 



1 '^l 



3 -f^y 



1 



(15) 



We define the (4x4x4) tensor r,f^\, „ and the len£rth-4 column vector r,^^ so that the 
above equation may be written in terms of Clebsch-Gordon coefficients of Sn , C^l^l'^lJ^ i 

22! ) and their coefficient tensors 



(see Section V.B of Ref. 

-Hi 

^J■2,^i 



^ _i (1) G - - I i_ (1) 



T, 



-Yd 



■M2 3 Vl,/^2,M2,/ ^M2 "I" 2 1 Vl^Ml "I" 1 Vl ■ 



(16) 



M2 



Therefore, the polynomial (//i)^^ may be written in the following compact form: 



'eff / ^ '^^^l,^J■2,^^3,^^ 

fll,fl2,l^3,T^ ^1,^2,^3 
At={0+,0-} 



(17) 



2. Derivation of the cubic A 



From Eqs. (jH]) and (I17p (in Sec. HID ip . we obtain the polynomial equation 



[A,i7o]$o= (i^O^ff^o. 



(18) 
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Solving this equation for the polynomial A, we obtain 



A'i>/i2,M3,''^5i,6,53 /i={o+,o-} 

where 

(1) A _ 3 ^l]_,^^2,^^3,k 

3 ' /il,/i2,At3,fc 



1 '^o± — ~, I 1 '^0± 



(20) 
(21) 



-t- ^ \^ 3 ^0±M/i 3 VO±M + 3 VA'Oit y I • V^^J 



Therefore, the first-order manybody wave function is obtained by multiplying the lowest- 
order wave function by A, a polynomial in q' given by Eqs f|T9l) and fl20l) : 

$i(q') = (l + 55A)$o(q')- (23) 

III. DERIVATION OF THE DENSITY PROFILE THROUGH FIRST ORDER 
A. Recap: The Lowest-Order Density Profile 



In Ref. l2l|, we derived the Jacobian- weighted density profile at lowest order by integrating 



over many of the degrees of freedom of the wave function, and transforming that integral 



from normal to internal coordinates. We arrived at (Eq. (68) of Ref. l2ll ) 



SiD)M,(r) = iV,/;^f exp -R r-|£ - VO f^ ] | . (24) 



The factor S{D) is the D-dimensional solid angle 24l |. 



2nf 

S(Z))^^. (25) 

and the quantity R (a number) is defined in Eq. (!6T!) below. 

Notice that the lowest-order Jacobian-weighted density profile is a Gaussian (normalized 
to A^) centered around r = k,{D) foo , the D oo configuration radius in oscillator units (see 
Eqs. (9) and (13) of Ref. [21]). The form of this Gaussian function is flexible in the sense that 
its shape depends on the values of two quantities, R and foo . However, this lowest-order 
Jacobian-weighted density profile is still limited to a symmetric shape about r = k{D) foo . 
The first-order corrections will add further flexibility by allowing for asymmetry. 
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B. First-Order Corrections 



The derivation of the first-order density profile is similar to that of the lowest-order 
density profile in that the same transformations are used to perform a change of coordinates. 
Integrals over the normal coordinates g'o+ , q'q- , [q\+\N-i , and [q'-j_-]]\f_i are transformed 

](JV-1) ) 

AT-l 

S = 



to , , S]^^ , and [Sl^ ^' ^^](n~i) , where 



-/ 



(26) 

1=1 

The first-order density profile is derived from the first-order wave function in a similar 
way to lowest order derivation. In Ref. 21 we showed that we can write 

S{D)N{r) 

■■■ 5/(r - r^) y<l>(y')]' H U^W,^^ ' (27) 

By simply substituting [g$i(q')]^ = (1+6^/^ Af [g^o{^)f for [g$(q')]^ in Eq. ([27]) we obtain 
the first-order density profile as 



OO POO 



X 



<5;(r - r^) (1 + + ^A^) [,<l>o(q')]' , (28) 



where 6f{r — Vi) is the Dirac delta function (differentiated from the inverse dimension, 6 , by 
the subscript /). Thus to first order in 5^/^ the Jacobian- weighted density profile is 



AT- /"OO /-OO '^M 

M(r) = Afoir) + ;^25V2 / . . . / JJ 

x5^(r-r;v)A [,$o(q')]' • (29) 



Substituting A from Eq. (|T9l) . we obtain 

AT poo poo '^M 

M(r)=Aro(r) + — -25^2 / .../ Yl U^W.h^A^- 

^[^) J~oo -/-OO ± ± 



X 



fll,^J.2,f^3,Tl 5l,52,53 /i={0+,0-} 



X [.$o(q')]' , (30) 

where /ii , /i2 , and /X3 each run over 0+ , 0^ , 1^ 1" , and 2 , and 1 < .^i , ^2 , (,3 ^ dfi 
where ci^ = 1 , — 1 and A^(A^ — 3)/2 for /i = 0^ , 1^ , and 2 respectively. When fi = 2 , 
I <n<II; otherwise 7^ = / . 



1. Normal Coordinate Integrals 

We simphfy Eq. f l30|) by defining the {P x P x P) tensor ^Mi^-^^iy^^i,,^ of integrals and the 
length P column vector iM,^ of integrals as special cases of the rank n tensor of integrals 



/oo /»oo ^ 
■■■/ n Wd[%]ii^fi'^ -'^N)qu,(l'u2---<lu„ [g'^oi^)? ■ (31) 
00 J —CO -« -I- n <- 1 



V=0±,l±,2 5=1 



Indexing 3M and \M by {fii,^i) rather than Ui (See Ref. |25|, Table III), we can write the 
density profile as a tensor contraction 

Afi{r) = Afo{r)+ (32) 



SiD) 



.^J.l,fl2,^J.3,'R- 5i.52,53 a'={o+,o-} 



Each element of the M tensor is an integral. 

Due to the presence of the Dirac delta function in 3M and iM we must treat integrals 
over normal coordinates involving differently from those that don't involve tn . As in 



Ref. 



21 



which derives the lowest-order density profile A/'o(r) , we note that appears in 
only four normal coordinates, i.e. in q'Q+ , q'^- , and (see Eqs. and 

(jl])). The remaining {N{N + l)/2 — 4) integrals, the overwhelming majority of the integrals, 
are easily evaluated using the integral identity 

^ 2./^uJ V 2 y 




which is zero when m is odd. Thus many terms are zeroed out; only those terms involving 
even powers of the normal coordinates which are independent of r^v contribute. 
Thus from Eqs. (l32l) and (l33l) we find 
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(34) 

where all primed indices range over 0^ and 1^, S^i o± equals one when /i' = 0^ but is 
zero otherwise, /i ranges over 1^ and 2 , and 1 < ^ < di± — 1 = N — 2 when /2 = 1^ , or 
1 < f < ^2 = A^(A^ - 3) /2 when fi = 2 . 



2. Clehsch- Gordon tensor contractions 



The relevant Clebsch- Gordon elements and sums in the above equation are 

r^O 0./ _ 1 

'-'111' — 

r^i ./ _ ,^0 1 ,/ _ ,^0 1 ,/ _ n 

^N-l 11' — N-l 1 — '-"1 1 N-l — U 

1 ,/ _ r^i 1 ,/ _ /^o 1 1 ,/ _ 1 (35) 

"-"Af-l N-l 1 — "-"Af-l 1 Af-1 — "-"l N-l N-l — ^ ^ ' 



^111,/ 

•-^Af-l Af-1 Af-1 



-(iV-2) 



^/N{N - 1) 



Af-2 



5^C||?'^ = iV-2 (36) 



Af(Af-3)/2 



^ (37) 



«=1 

^-2 ^ _ 2 

g^{f"-''= Viv(iv-i) 

Af(Af-3)/2 

E C7|-_,'^ = 0. (39) 

Equation fl39|) simply follows from the fact that since the 2 indices are saturated, there is 
nothing to couple the 1 index with to form a scalar irrep. . 

Equation ( !34l) may be further simplified by defining two tensors which hold the coefficients 
of the elements of ^M^^^^f'l ^ ^ and iM^ ^ tensors. We define a 4 x 4 x 4 tensor, '^E^' y , 

to be the non-zero coefficients of the elements of sM^^^'^f^ ^ ^ (which are cubic in tn) 
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We emphasize that the above equation is a simple elemental multiplication, no summation 
is implied. Using the above Clebsch-Gordon elements, we obtain 



^i?o±,o±,o± - T^,o±,o± (41) 

^-E'i±,o±,o± = ^-So±,i±,o± = ^-^'Oit.Oit.iit = (42) 

^£^i±,i±,o± = ^i^,i±,o± (43) 

'i5i±,o±,i± = ^I'i.oi.ii (44) 

^^o±,i±,i± = ^o±,i±,i± (45) 

In Eqs. (HT]) - (H^ . each ± associated with a sector fi is taken to be independent of the ± 
associated with the other two sectors. 

We also define a length- four column vector ^i?^/ to be the non-zero coefficients of the 
elements of iM^^ (which are linear in r^)- 

^^M' = X] 2^ ( ^3 + ^3 ^a!/^',m,7 + ^]i'^^,tl,tl,I ) '^UX, + '^^^''0* ^l^^M' ' (4"^) 



I.e. 



(iv-2; 




2 






2) 1 


2 




iV(iV 


-3) 1 



ip _ l^v - ^) i ^ A , A , A ^ 

-^Oit — ^ ^ ^ V l+l+0± "1" ^1+0±1+ ^0±1+1+;' 

/A I I \ 

Vl-l-0± "I" ^1-0±1- ^Oil-l-J 

A , ('^220± + ^20±2 + ^0±22) + ^^^0± (48) 

Ip _ 1 /"W^A , (1) A , (1) A 

- 2^ iV(iV-l) ^l+ V 3^1+1+1± + 3^1+ + 3^1±1+1- 

I ^~ ^ 1 Al)^A , (1) A , (1) A \ /.qN 

+ 2v/ iV(iV-l) ^ I ^ "^^-^-^^ + ^ + ^ ^ • ^^^^ 



Using Eqs. fHOl) and (l47j) . Eq. flMl) may be written in a simpler form as 
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3. Transformation of the integrals to symmetry coordinates 



The elements of the M tensors are integrals over the normal coordinates. We use the T 



transformation of Ref. 



21 



where 



q' 



T a' 



(51) 



to write these M tensors in terms of integrals over the coordinates of the four-element vector 



a 



' N 
b' 



(52) 



where 



r's 

aim 



(53) 



I.e. 



^ d, d,. 



Ml Mo M 



(54) 



ll,l2,. ..,ln 



where T is given in Eqs. (53), (54), and (55) of Ref. l2l|. Applying the analysis of Section VIII 
of Ref. [21 to the ri-M.ixi2--in tensor, we arrive at 



Jt\/uJq+ UJq- UJ-1+ UJi- 



«2 ■■■ ^n 



00 /•OO /"OO />oo 



00 «/ —00 J —00 J —00 



5f{r - fli) 



(55) 



tti-^ai^ ■ ■ ■ tti^ exp (— -K'o af — 2aiK'^b' — b'^Kb'^ dai Sb' . 



where Kq , K , and K are defined in Eqs. (59), (60) and (61) respectively of Ref. 



21 



and Jt 



is the Jacobian of the transformation from q'Q+ , q'^ , [^'1+]^^+ and [q\-]d^- to the internal 
coordinates a' , Eq. (55) of Ref. 



21 

Substituting Eq. fl5^ into Eq. fl50|) we obtain the density profile in terms of integrals over 
internal coordinates, i.e. 



M(r) =Afoir) + 



2(5^^ N 



il,«2,«3 « / 



(56) 
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where the (length 4) column vector iSj and the (4x4x4) tensor sSjj-^A; are 

iSi = ^ ^i?^' T^i^i (57) 



4- Evaluation of the integrals 

To perform the integrals of Eq. fl55|) for sA^jiijjs ^'^d lA^j in Eq. fl56l) . we define the 
four-component vector 

Eq. (155|) can be written as a series of derivatives of a term proportional to A/o(r) 



V det K R 

nA<ili2...i„ = , (60) 

/-i a \ /-I a \ /-i a \ /■°° f"" r , 



exp (-Ko al - 2aiK^b' - b'^Kb') dai cfb' . 



where we have also used Eq. (67) of Ref. l2l|, 



R = ^0-^2"^^.- 4 ^ _ KTj^-iK^ . (61) 
Upon using the integral identity 

/OO POO 
■■■ exp {-b'^Kb' - 2ai K^b') (Tb' 
-OO J — OO 

exp [al K^K-^K) , (62) 



-CXD J — OO 

TT 2 



V det /C 
in Eq. (!60l) we obtain 



(63) 



X 
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(64) 

which yields 

M l^i5 l^ie X ■ ■ ■ X l^in ^' (65) 

_^ ^ 2^n i2 2^i3 i4 X ■ ■ ■ X 2^i„_i in ) when is even , 

2Xii i2 2Xi3 X ■ ■ ■ X 2Xi„_2 iXi„f' ) when n is odd , 

where the Xs are elements of 

= -^Vv(-i?)= f \ I (66) 

1 1 [ \ 

2X = --Vv®iX = - , (67) 

2 2 ^0 y 

iX is four-dimensional vector and 2X is a 4 x 4-dimensional matrix. The C operator acts 
on each term in Eq. fl^^ to produce a sum of terms over all the distinct combinations of 
indices. 

From Eq. ( l65i) we can now evaluate each element of the Ai tensors: 



3-^iii2i3 — \l c I r^\2 ( l^ii 1^42 l-^is '"'^ + ( 2^n 42 l^ia + 2^ii 13 1^12 + 2^12 13 I'^^ii ) • 



5k{D) 



TT 



(69) 



C. Result: first-order density profile 

1. density profile in displacement coordinates 

Using Eqs. fl68l) and fl69l) in Eq. fl56|) yields the density profile 



M(r) = -^^^J^^i^ + '^^(^1^" + ^3r'^))) exp(-i?f'^) , (70) 

15 



where 



f'(r; D) = VD (^-^ - foo) . (71) 



The coefficients Ai and of the polynomial are 



Ai — 2 ( ^ ^ iXj + ^ ^ ( 2^ii 42 1-^13 + 2^11 13 1^12 + 2^12 is I'^ii) 12 is 1 

\ i «i,*2,j3 / 



(72) 



A3 — 2 i^ii 1^12 1^13 S'^'ii 12 «3 ' (73) 

where 

iS, = J] T,, ^i?, (74) 

At={0±,l±} 

3^i,j,k = ^ ] ^ ] ^ ^ ^/il,i^/X2,j^/i3,fc -^Mi,/i2,M3 5 C^^) 

/il = {0±,l±} M2={0±,1±} /i3={0±,l±} 

iX,= <^ ^ , (76) 



and 



-(/C-iK)i_i 1< i < 4 

ix,efc_ie,_i , (77) 



where 



(78) 



©m = when m < 
= 1 when m > 1 . 

The ffist-order density profile has the form of a cubic polynomial multiplied by the lowest- 
order density profile. Note that the density proffie is a function of the coordinate r, which 
is not the dimensionally-scaled internal displacement coordinate f'{r]D). Thus one must 
make the following substitution to obtain the density proffie as an explicit function of r: 

f'(r;D) = r^(^-foo) . (79) 

We have derived the A^-body density profile through ffist-order in 5^/^ . This density 
profile, Afi (r) , of Eq. ( ITOl) includes the full interactions of particles through ffist order 
in the perturbation series, exactly and analytically. The details of the interactions at each 
order are folded into R, Ai , A3 , and foo . We note though, that we have neglected the 
next-order term (order 6) in Eq. (!28l) . 
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2. Density "profile in oscillator units 

For the case of a system under harmonic confinement (such as the present trapped Hooke's 
law gas or a BEC in a parabohc trap) we may choose a oscillator-unit scaling k,{D) = D'^aho, 
where 

1 

O/io — j^^j2 ^ho ■ (80) 

Therefore 

K{D) = y/Daho (81) 



f = (— -v^foo), (82) 

O-ho 



and defining 



r 

'"osc 



a-ho 

in oscillator units, we obtain the Jacobian-weighted density per particle 



^exp (^-R (tosc - VDt^^ ^ 
M (rose) X = (^1 + 53 (^Ai (^rosc -VDfJj+ (^r^sc - Vd 



3 

r, 



(83) 



ATo (rose) X = \l- exp ( -R ( rose -VDr^) ] (84) 



y<^fo{rosc)y<^■ (85) 



IV. TEST APPLICATION: HARMONICALLY INTERACTING PARTICLES UN- 
DER HARMONIC CONFINEMENT 

In this paper we have derived the Jacobian-weighted density profile through first order 
from the wave function through first order as an example of an observable which may be 
derived from the A^-body, interacting wave function. We test this general formalism for the 
density profile, by comparing it to the density profile of the analytically solvable system of 
N , harmonically-interacting particles in a harmonic confining potential with Hamiltonian 



The exact analytic density profile for this system is is independently derived in Appendix [B] 
(see Eqs. (IB48P and (IB54p ). and this is expanded through first order in 5^/^ to yield the 
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exact density profile through first order (see Eq. (1B70P ). This analysis shows that the 
density profile for any or interaction strength follows a universal curve when a simple 
scahng is applied to the radial variable (it should be noted that this is not true of the wave 
function) and is given by 

A/'(^eff) = ^ f^s' exp {-D fls) , (87) 

where 

^eflf = VXff r , (88) 



and 



Aeff= A= ^l + iVA2, Xp = Up/ujt (89) 



AfiK,) = ^^fir'). (90) 



Expanded to first order Eq. (187|) gives 

mK,) = (l + 5^ V2 - K,) ) (^) ' exp (-2 rt,) . (91) 

This scaled density profile for D = 3 (5 = 1/3) is plotted in Fig. [H One readily sees the 
improvement obtained at first order, confirming the efficacy of this approach to the general 
confined A^-body problem, which may be systematically improved by going to higher orders. 

The general theory developed in this paper for the density profile involves no such har- 
monic interaction specific scaling since it's applicable to any interparticle potential, not just 
harmonic interparticle potentials. Consequently in Figs. [2] and [3] we plot the density profile 
for D = 3 {6 = 1/3) without this harmonic-interaction specific scaling for two very different 
interparticle interaction strengths. Both are for = 10, 000 particles, but Fig. [2] features 
strongly attractive interactions, while Fig. [3] is for a repulsive interaction just below the 
dissociation limit. In the former case the system is tightly bound and very compact. In 
the latter case the confining potential is barely able to hold the system together against the 
combined effect of the repulsive interactions, and the system is very extended. 

The density profile derived from the general A^-body formalism developed in this paper. 



and implemented in Mathematical:^ code 



27( 1 , is indistinguishable from the density profiles 



derived from the exact independent solution of the harmonically-interacting system. The 
agreement between the general formalism, which uses a perturbation series invariant under 
Sn , and the direct density profiles of Eqs. (IHTl) and (iQljl obtained in Appendix [Bl confirms 
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the correctness of the general formahsm developed in this paper, and its implementation in 
MathematicaU& code. 




V. SUMMARY AND CONCLUSIONS 

While the resources required to solve classical systems scale as a polynomial with the 
number of particles, N , allowing calculations to be performed involving large numbers of 
particles, the situation regarding large- iV quantum systems is not so felicitous. In this case 
the resources required scale exponentially with , making calculations for large- quantum 
systems a far more formidable challenge, unless N and the interparticle interaction strengths 
allow an approximation, such as the mean-field approximation. 

In a series of papers, we have been developing an approach to the general interacting 
quantum A^-body problem, which while essentially analytic, makes no assumtions regarding 
the form or strength of the interparticle interactions. This approach derives the interacting 
A^-body wave function, from which any observable quantity can in principle be derived. In 
the process, collective, normal mode coordinates are derived revealing the nature of the 
microscopic motions of the particles of the system for any interaction. Using this wave 
function, properties such as energies and density profiles have been derived. 

The method involves expanding the system in inverse powers of the spatial dimension 
D . At large D , systems exhibit a point group structure of a far higher degree of symmetry 
(isomorphic to S^) than is possible in three dimensions, allowing group theory and graphical 
techniques to be used to tame the exponential A^ scaling, leading to an essentially analytic 
solution at lowest order in the wave function. 

More recently, in a major development of this approach, this method has been succesfuUy 
extended, analytically and exactly, to first-order in the wave function, and in principle the 
techniques developed to do this can be extended to yet higher orders. 

This paper is the first application of this first-order wave function to the derivation of 
physical property of interacting A^-body quantum systems, namely the density profile to 
first order. In a test of this theory, the derived first-order density profile is tested on an 
exactly solvable model, namely a system of harmonically interacting particles in a harmonic 
confining potential. The harmonic interparticle interactions may be attractive, or repulsive, 
and if the interactions are sufficiently repulsive the system will dissociate despite the presence 
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of the harmonic confining potentiaL 

The general theory developed in this and prior papers, agrees with the exact results for 
this system obtained from the independent solution, showing strong convergence to the exact, 
three-dimensional result for both strongly attractive interactions and repulsive interactions 
just below the dissociation threshold. 



While this paper (as well as Ref . llTl ) has focused on the harmonically-interacting system in 
a harmonic confining potential, the theory is not limited to these systems and in past papers 
we have examined other systems at lowest order, such as the Bose-Einsten condensate, and 
quantum dots. While the lowest-order approximation for the BEC captures the behavior of 
the system for a range of and interaction strengths, when or the interaction strength 
is larger than this range, the lowest-order density profile increasingly does not have the 
fiexiblity needed to capture the behavior of the system. The first-order result does not have 
this limitation and so it is very desirable to apply the first-order density profile derived in 
this paper to strongly interacting BECs as well as other strongly interacting systems. Of 
particular note is the fact that the functional form of the first-order density profile of Eq. (l85l) 
admits structure (wiggles) indicating the onset of crystallization/fermionization. Although 
such behavior is not seen for the long-range harmonic interactions examined in this paper, 
other systems do exhibit such transitions. 

While we focused in this paper on density profiles, the theory derives the exact, first- 
order wave function from which any observable may be derived. Also excited states may be 
derived, and the theory needed to extend these results to higher-angular-momentum states 



has been set up [28 



29| 
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APPENDIX A: BINARY INVARIANTS 



In Ref. 
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we introduce binary invariants which are binary tensors invariant under the 



Sn particle interchange group. 



1. Introducing Graphs 

Definition 1 A graph Q = {V,E) is a set of vertices V and edges E. Each edge has one 
or two associated vertices, which are called its endpoints. 

For example, J is a graph Q with three vertices (or "dots") and two edges (or lines). 



We allow our graphs to include loops and multiple edges jSOj. A graph contains information 
regarding the connectivity of edges and vertices only: the orientation of edges and vertices 
is insignificant. 

We introduce a mapping which associates each tensor element with a graph as follows: 

1. draw a labeled vertex (Q i) for each distinct index in the set of indices of the element 

2. draw an edge ( LJ) for each double index (ij) 

3. draw a "loop" edge (O for each distinct single index i 

For example, the graph corresponding to the tensor element ^2'*Qi|ij) under this mapping is 

Two graphs with the same number of vertices and edges that are connected the same 
way are called isomorphic. The action Sn group interchanging particle labels, only connects 
tensor elements with isomorphic graphs; two elements with heteromorphic graphs are never 
connected by the Sn group. Now consider a tensor, for which all of the elements labeled by 
a a single isomorphic set of graphs are equal to unity, while all of the other elements labeled 
by graphs heteromorphic to this single set of isomorphic graphs are equal to zero. We we 
term this tensor a binary invariant, [B{Q)]^-^^u^^,,, since it is invariant under the Sn group, 
and we label it by the graph Q; sans particle labels at the vertices, for the non-zero elements 
all of which are equal to unity. 

Denoting the set of unlabeled graphs for each block as GxiX2...x„, where n is the rank of 
the tensor block (and therefore the number of edges in each graph in the set) and X is r or 
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7 , we have 



Grr = {OO, 0} 

Gyr = {0-,0 (Al) 

= {o, 



Gr = {0} 

G, = {-} 

Grrr = {Q^, gO,000} (A2) 

= {8^,cx3,o-0,oo,oj:)} 
G777 = {e, ZI'O-' J-'U>20'S} 

Each of these graphs denotes a binary invariant. 



APPENDIX B: DENSITY PROFILE OF HARMONICALLY INTERACTING 
PARTICLES SUBJECT TO A HARMONIC CONFINING POTENTIAL: AN IN- 
DEPENDENT SOLUTION 



The Hamiltonian of the harmonically-interacting model system of identical particles is 

N r 02 



-4 E 



drf 



i<j 



(Bl) 



1. The exact AT-body Wave function 



Making the orthogonal transformation to center-of-mass and Jacobi coordinates 
R=^p^ and p. = -^gr,-ir,„^ 



where 1 < i < N — 1 , the Hamiltonian becomes 



2 



92 



\ 1 



i=l 



92 



2 2 



the sum of A*" , dimensional, harmonic-oscillator Hamiltonians, where 



(B2) 



(B3) 



(B4) 
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The ground-state solution of the Schrodinger equation 

H^ = E'^ (B5) 
is the product of harmonic-oscillator wave functions 

^{R, {pi}; D) = ij{R; Ut, D) J] ^(p,; cui^t, D) , (B6) 

i=l 

where ^/^(r; uj, D) is the D-dimensional, harmonic-oscillator, ground-state wave function 



_D 



i,{r- u, D) = \l-^ exp [-^ ) (^7) 
and iplr; uj, D) satisfies the normalization condition 

{r; uj, D)fr^^W^r = 1 (B8) 



(B9) 



so that 

7V-1 



/oo J- 
■~ i=i 



2. The exact A'^-body density profile 

Since the wave function is a completely symmetric function under any permutation of 
the particles, one can write the Jacobian-weighted density profile as 

N-l 



Af{r) 



/oo ^ ' J- 

[^{R, {p,}; D)]^ 5(r - r^) J] rf^p, d^R , (BIO) 
i=i 



where from Eqs (1B6P and (]B7p . the unweighted wave function is 



JV-l 

D 



2 N-l 



{P&, D) = f-^ exp i-^-lR-) (^^j n -P • 

Since the transformation from single-particle coordinates, ri , . . . , r^v , to Jacobi/center-of- 
mass coordinates, {pi}/R, is orthogonal. 



/ D \ 












II 




\-n2 i 







^{r) =1^1 1^1 M(r) = r^-V(r) (B12) 
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where p(r) is the unweighted density profile, and 



M(r) = J d^n d^r2 . . . d^rNS{r-rN) exp {-{cut - uJi^t)R^) exp (-cui^t + I] P: 



(B13) 



Defining 



jvr = ri © ra 



r2 



N 



J = ® Id , 



(B14) 



(B15) 



where Id is the D x D dimensional unit matrix in the D-dimensional coordinate space and 
Jiv is the N X N dimensional matrix of elements equal to one in the particle label space, 
then we can write 

R^^^nV^nJnT (B16) 

7V-1 TV 



1=1 



i=j 



where 



nI = In® Id 



(B18) 



and Iat is the N x N dimensional unit matrix in the particle label space. Since J d^VN — 
J r^~^drN J dfijv, where J dfijv is the integral over the D-dimensional solid angle, we can 
write 

M{r) = J d^ri d^r2 . . . d^VN S{r - tn) exp ^_^t^L._p^ nv'^nJ nv) exp {-cui^t nv'^nI nv) 



S{D) r^-^ exp 



{Ut + (iV - l)tU^nt) 

N 



X 



J d^ri d^r2 . . . d^r^-i exp (-Ar_ir^A ^-ir - 2 N-ir) , 



(B19) 



where 



S{D) 



„ D 

2%- 



r(f) 



(B20) 



is the £)-dimensional solid angle. 



N-ir = ri © r2 



rjv-2 ® r-jv-i , 



(B21) 
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A = cuint N-iI + ^ * ^ '^^^ N-lJ (B22) 
B = ^J^^l^pA 4 r (B23) 



J = Jw_i ® /d , (B24) 

Jat-i is the (A^ — 1) x (A^ — l)-dimensional matrix of elements equal to one in the particle 
label space of the first — 1 particles, 

N-iI = In~i ® Id , (B25) 

and Iat-i is the (A^ — 1) x [N — l)-dimensional unit matrix in the particle label space of the 
first A^ — 1 particles, 

el = le ® , (B26) 

where Ic is the (A^ — l)-dimensional column vector of elements equal to one in the particle 
label space of the first A^ — 1 particles so that 

N-l 

cl^N-ir = (B27) 

k=l 

Using the result 

/OO -2; 
c/"6exp(-6^A6-2B^6) exp (B^A~^ B) (B28) 

OO VdetA 

in Eq. (]B19p . with the identification n = D{N — 1) and b = r , we obtain 

M(r) = S{D) r^-i exp ( J^t + (N - l)uj,,,,) vr^^ ^^^g) 

V ^ /V det A 

Thus we need to evaluate det A and A^^ . First A^^ . One has that {U <^V)^^ = U^^^V^^ . 
Thus writing 

U = a In-1 + P Sn-1 , (B30) 
where a = cuint and /3 = ^ ^^^^ using the closed algebra 

Jjv-i In-1 = (N-l) 
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we obtain 



= (« I;v-i + /? Jtv-i)"' = - (in-1 , Jtv-i ) • (B32) 

a \ a + [N — l)p / 

So with V = Id , one finds that 

A-' = — ( In-. - J^-i) ® Id • (B33) 

^^int V {N -l)UJt + UJint J 

To evaluate det A we note that 

det {U®V) = (det C/)'^^ (det V)'^^ , (B34) 

where d^/ and dy are the dimensionahties of matrices U and V respectively. The determinant 
of Id is simple enough and equals unity. For det U of Eq. (1B30P , we show that 

det U = det [a I^v-i + /? J^v-i] = a^^^^^ (a + (A^ - 1) /3) (B35) 

Proof: 
We have 

det [a In-1 + P Jjv-i] = ||^^|| , (B36) 



where E is the diagonal matrix of eigenvalues of V, and \\E\\ is the norm of E , the product 
of its eigenvalues. Now sin 
eigenvector of U satisfying 



of its eigenvalues. Now since Jat-i = Icl^ , x^^i = Ic is seen to be a normalized 



U xn-i = En-i xn-i , (B37) 

where 

EN^i = a + {N -1)(3. (B38) 

The remaining N — 2 eigenvectors, Xj , \f 1 < J < iV — 2 , are orthogonal to xn-i so 
In-1 Xj = ^ ^ 1 < j < ^ — 2 , from which it follows that 

U Xj = Ej Xj , (B39) 

where 

Ej = a V 1 < J < - 2 . (B40) 
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Putting it all together we arrive at Eq. (lB35p . Q.E.D.. 
Using Eq. flB35|) in Eq. flB34|) we arrive at 



det A = (a + (AT - 1) (B41) 

D(N~2) 

- ((iV-lV, + a;i„0^ (B42) 



As for exp {B^A-^ B) , Eqs. flB23l) and flB33|) . along with 



li Ijv-i Ic = (iV - 1) (B43) 

iJJiv-ilc = (iV-1)', (B44) 

give us 

exp [B-A-^ B) = exp f ~ ^^^1""^' ^'1 ' ^^45) 
Thus from Eqs. flB29D . fiBiB . and fiBiSD we obtain 

TT 2 /V 2 

M(r) = ^(D) r^-i = exp {-uj,^ r^) , (B46) 



where 



Finally using Eqs. (1B20D . flB46D . and (IB47D in Eq. (IBT2D we obtain 



D 
2 



Ar(r) = ^^^^^ r^-^ exp (-A^ff^, r^) , (B48) 



where 



Note that 



(B49) 



L^int + (AT - l)cUt 2f^ 

N{r)dr = 1 . (B50) 



3. Dimensional Expansion of the Number Density 

a. Dimensional scaling 

To obtain the dimensional expansion of the density profile of Eq. ( ]B48|) . we first need to 
regularize the large-dimension limit by dimensionally scaling the parameters and variables. 
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As in Ref. Il9|, we define dimensionally scaled frequency uot and tlie dimensionally-scaled, 
oscillator-scaled radial variable f, 




from which we derive 



(B51) 
(B52) 

(B53) 
(B54) 

Af{f)df = 1 . (B55) 

Equation (1B54|) implies (as does Eq. (1B48|) ) that up to a scaling, the density profiles for 
harmonically-interacting particles in a harmonic-confining potential follow a universal curve 
for any or interparticle interaction strength in oscillator units Xp = uOp/ ujt , where 

(B56) 



Thus the number density in dimensionally-scaled coordinates is 

m = ^-^j^ f^-' exp (-Ae. D f^) , 

where 



A 



A^A 



l + ATA^ and Aeg = ^ ^ _ ^ 



This is simply seen by scaling the dimensionally-scaled radius r : 

fes = v^A^ r 

and scaling the wave function by the multiplier 1 / V%ff from the change of variables 

1 



dr 



dr, 



eflf 



(B57) 



(B58) 



cfT 



which gives the equation for the universal curve of the density profile as 

-^(^^eff) = ^ fff' exp {-D rig) , (B59) 

where 

POO 

/ ^f{f,s)dns = 1 . (B60) 

As -D — s> oo , the number density of Eq. (lB54p becomes more and more strongly peaked 
at (feff)oo determined from 

'^^^''^'^ 0, (B61) 



eff 



D^oo 



and thus the peak of the number density occurs precisely at the large-D radius parameter 

' . _ 1 



(B62) 
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b. Series expansion 



As in Ref. [l9|, we introduce the dimensionally-scaled displacement coordinate 



so that 



Thus 



where 



Tea = (r'efr)oo + S'^Ks = Tjir'es + D^r^s) 

D2 



D2 



/oo 

To derive the dimensional expansion of Eq. (]B65p . let's first consider expanding 
DUfcs) oo)^"^ ■ We derive 



D-l 



X exp 7—^ exp (- r^gj . 



00 



Likewise we also have 

exp (- (f^ff + I^^(fcfr)oo)') = exp exp (^-d'^^^^^ exp (-f^^) . 

Equations (IB62D . (iBeTj) and (|B68|) . along with 



2Vexp(^: 



0(5) 



r (f ) ^/^^ 

give us the result we are after, namely 

M(r'^^) = (1 + 5^ V2 - rl,) + 0(5)) 0) ' exp (-2r- ) 

where through order ^2 the normalization condition 

/oo 
A/'(r:ff) rfr-f, = 1 
-00 

is still satisfied. 
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As we noted after Eq. (IBSSp . the density profile for harmonically-interacting particles in a 
harmonic confining potential follows a universal curve for any or interparticle interaction 
strength Xp = Up/uf. Although the density profile has this property, the same cannot be 
said for the wave function for A^ , harmonically-interacting particles in a harmonic confining 
potential. For example, many terms in the wave function through first order in 6^ (see 
Eqs. (24), (25) and (26) of Ref. I22I) are zero for the free trap (Ap = ^ A = 1 Aes = 1 , 
^00 = l/"\/2, and 700 = 0). Thus there is no simple scaling between the wave function 
for non-interacting particles in a harmonic confining potential and the wave function for 
harmonically-interacting particles in a harmonic confining potential. 
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FIG. 1: Scaled density profile at D = 3 for N harmonically-interacting particles under harmonic 
confinement in oscillator units of the confining potential. The short dash curve is the lowest-order 
density profile, while the solid curve is the density profile through first order. The long-dash curve 
is the exact result. The scaling factor, \/XeS is explained in Appendix [Bl 
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FIG. 2: Unsealed density profile at D = 3 for N = 10,000 particles under harmonic confinement 
with strong attractive harmonic interactions (Ap = 100) in oscillator units of the confining potential. 
The short dash curve is the lowest-order density profile, while the solid curve is the density profile 
through first order. The long-dash curve is the exact result. The parameter Ap , as explained in 
Appendix [B1 is the interaction frequency squared in oscillator units of the confining potential. 
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FIG. 3: Unsealed density profile at D = 3 for N = 10,000 particles under harmonic confinement 
with repulsive harmonic interactions (Ap = —1/10,000 + lO"^'') in oscillator units of the confining 
potential. The system is just below the dissociation threshold and very extended. The short dash 
curve is the lowest-order density profile, while the solid curve is the density profile through first 
order. The long-dash curve is the exact result. The parameter , as explained in Appendix [Bl is 
the interaction frequency squared in oscillator units of the confining potential. 
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